load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

  varname = "NUMICE"

  plotname = "profile_"+varname+"_acpd_v3" 
  plotform = "eps" 

  ;;.......................
  ;; read data 
  ;;.......................

  fna = "../diag_cubic_fixq.nc"
  fnb = "../diag_milled_fixq.nc"
  fna = "../diag_cubic_acpd_v3.nc"
  fnb = "../diag_milled_acpd_v3.nc"
  fnc = "../sgp_20100401_shaocheng.nc"


  fla = addfile(fna,"r")
  flb = addfile(fnb,"r")

  xa = fla->$varname$
  xb = flb->$varname$
  xt = fla->T
  xp = fla->PS
  hyam = fla->hyam
  hybm = fla->hybm

  nt   = dimsizes(xa(:,0,0,0))
  nlev = dimsizes(xb(0,:,0,0))

  v1 = xa(:,:,0,0)
  v2 = xb(:,:,0,0)
  vt = xt(:,:,0,0)
  vpp = xp(:,0,0)

  rho = vt
  pre = vt

  do ik = 0,nlev-1
     pre(:,ik) = (/ doubletofloat(hyam(ik)) + doubletofloat(hybm(ik)) * vpp(:) /)
  end do

  r = 287.042 ;; J/kg/K.

  rho = pre / (r * vt)

  ;;.......................
  ;; #/kg -> #/L 
  ;;.......................

  v1 = v1 * rho / 1.e3
  v2 = v2 * rho / 1.e3

  ;;.......................
  ;; define new vars 
  ;;.......................

  vv = new((/2,nlev/),"float") 

  printVarSummary(v1) 

  ;;v1 = where(v1.gt.1.e-9,v1,-999) 
  ;;v1@_FillValue = -999 
  ;; 
  ;;v2 = where(v2.gt.1.e-9,v2,-999) 
  ;;v2@_FillValue = -999 

  vv(0,:) = dim_avg(v1(lev|:,time|:)) 
  vv(1,:) = dim_avg(v2(lev|:,time|:)) 

  vv = where(vv.gt.1.e-9,vv,-999) 
  vv@_FillValue = -999 

  ;;.......................
  ;; plotting
  ;;.......................

  wks   = gsn_open_wks (plotform,plotname)                   ; open workstation
  
  res                   = True                       ; plot mods desired
  res@tiMainString      = "" 
  res@trYReverse        = True                       ; reverse Y-axis
  res@xyDashPatterns    = (/ 0, 15 /) 
  res@xyLineColors      = (/ "red", "blue" /) 
  res@xyLineThicknesses = (/ 4.0, 4.0/) 
 
  res@pmLegendDisplayMode    = "Always"              ; turn on legend
  res@pmLegendSide           = "Top"                 ; Change location of 
  res@pmLegendParallelPosF   = .70                   ; move units right
  res@pmLegendOrthogonalPosF = -0.25                  ; more neg = down
  res@pmLegendWidthF         = 0.12                  ; Change width and
  res@pmLegendHeightF        = 0.10                  ; height of legend.
  res@lgLabelFontHeightF     = .02                   ; change font height
  res@lgPerimOn              = False                 ; no box around
  res@xyExplicitLegendLabels = (/" Cubic"," Milled"/) 

  res@mpShapeMode  = "FreeAspect"
  res@vpWidthF      = 0.4
  res@vpHeightF     = 0.6

  res@tiXAxisString   = "NUMICE (#/L)" 
  res@tiYAxisString   = "Pressure (hPa)" 
print("vv : "+vv) 
 
  res@trXLog = True
  res@trXMinF = 0.01
  res@trXMaxF = 10.

  plot  = gsn_csm_xy (wks,vv,v1&lev,res) ; create plot

end



